Quantum Monte Carlo calculation of the finite temperature Mott-Hubbard transition 
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We present clear numerical evidence for the coexistence of metallic and insulating dynamical mean 
field tlieory(DMFT) solutions in a half-filled single-band Hubbard model with bare semicircular 
density of states at finite temperatures. Quantum Monte Carlo(QMC) method is used to solve the 
DMFT equations. We discuss important technical aspects of the DMFT-QMC which need to be 
taken into account in order to obtain the reliable results near the coexistence region. Among them 
are the critical slowing down of the iterative solutions near phase boundaries, the convergence criteria 
for the DMFT iterations, the interpolation of the discretized Green's function and the reduction 
of QMC statistical and systematic errors. Comparison of our results with those of other numerical 
methods is presented in a phase diagram. 



PACS numbers: 71.30.-t-h, 71.27.+a, 71.28+d 



The Mott-Hubbard Transition [yg] is an outstanding 
open problem in condensed matter physics. The theo- 
retical progress in this field has been impeded by non- 
perturbative nature of the strongly correlated system. 

In recent years considerable progress has been made 
with a single band Hubbard model by using DMFT ap- 
proximation which becomes exact in the limit of infinite 
coordinations (z = oo) |3J. The most striking point 
of these works is a discontinuous metal-insulator tran- 
sition(MIT) at finite temperature T which takes place 
at the critical interaction strength line Uc{T). This line 
lies inside a coexistence region where two phases, one 
metallic and one insulating coexist. The coexistence re- 
gion is bounded by two lines Uci{T) and Uc2{T) and 
Uci{T) < U,{T) < UMT) §. 

While this picture is now generally accepted [H by all 
groups working on the problem, the origin of the con- 
troversial results over the coexistence region in a single- 
band Hubbard model at finite temperature within QMC 
method, and more generally the difficulties in carrying 
out the DMFT-QMC algorithm |||-0], have not been 
clarified. This is an important problem, since DMFT- 
QMC is currently being applied to many similar models 
with more complicated orbital and spin structure. 

In this paper we provide clear evidence for the coexis- 
tence of metallic and insulating solutions in a half-filled 
single-band Hubbard model with bare semicircular den- 
sity of states at finite temperature within DMFT-QMC 
calculation [g[. We discuss several aspects of DMFT- 
QMC method which have to be taken into account in 
order to get the reliable results near a coexistence re- 
gion. First, we identify the source of the main difficul- 
ties in running DMFT-QMC near the coexistence region: 
uncontrolled errors in DMFT-QMC in the region where 
critical slowing down is dominant. The critical slowing 
down of the DMFT iteration follows from the Landau 
analysis M and can only be avoided by choosing a differ- 
ent rearrangement of the DMFT equation such as the one 
suggested by the Newton method. In this paper, however 
we use this critical slowing down as an indicator which 



allows us to locate boundaries of the coexistence region, 
and to fully confirm , using QMC methods, the predic- 
tions of the Landau analysis. We discuss how to minimize 
errors in DMFT-QMC calculation by using correct in- 
terpolation of the discretized Green's functions in direct 
Fourier Transformation(_FT) as well as a proper treat- 
ment of discontinuities of Green's functions in inverse 
Fourier Transformation(iFT) . Also extrapolation of data 
to the limit of dT{— (3/ L) ^ is used to reduce Trotter 
errors. Finally we obtain reliable coexistent insulating 
and metallic solutions within DMFT-QMC method and 
estimate the location of the first order MIT line Uc{T). 

In the limit of infinite dimensions a single-band Hub- 
bard model is mapped onto an Anderson impurity model 
supplemented by a self-consistency condition H . DMFT 
equation of this model. 



t^G{iuj)[A] = A{toj), 



(1) 



is solved iteratively within QMC calculation, where 
A(zw) and G'(«a;)[A] are a hybridization function and a 
local Green's function of the impurity model. Our energy 
unit is set to be i = 1/2 which comes from a scaled hop- 
ping amplitude, t/y/z, in an infinite dimensional single- 
band Hubbard model. 

Within DMFT-QMC calculation, the coexistence re- 
gion should be explored with a great care due to the 
following two reasons: critical slowing down and errors 
in DMFT-QMC. The critical slowing down is the criti- 
cal phenomena which manifests itself by slowing down of 
the convergence rate of a numerical method. The critical 
slowing down at a phase boundary can be proved with 
the argument based on Landau free energy |^ : no conver- 
gence at an inflection point of Landau free energy within 
iteration method. At a boundary of the coexistence 
region one solution never converges while another one 
converges fast. However around Mott endpoint (where 
Uci and Uc2 merge) both solutions converge slowly be- 
cause the narrow coexistence region around Mott end- 
point brings two boundaries close enough so as to make 



a whole coexistence region numerically unstable as was 
pointed out in ref . Q . Unless carefully taken care of, fluc- 
tuation of a solution originating from errors in DMFT- 
QMC may result in undesirable transition of one solution 
to another one in DMFT iteration around Mott endpoint 
where critical slowing down is dominant. 

In this paper we show that one can use the existence of 
critical slowing down as an indicator of boundaries of the 
coexistence region. Near a coexistence boundary, Uci{T) 
or Uc2iT), metallic and insulating solutions show differ- 
ent convergence patterns within DMFT iterations: one 
converges very slowly while the other converges in a few 
iterations. A convergence criterion is applied to deter- 
mine required number of iterations to get a convergent 
solution. We locate the coexistence boundaries at which 
such iteration numbers diverge. 

The convergence of a solution is determined by mon- 
itoring of the most fluctuating component ImTiiiirT), 
imaginary part of self energy at the first Matsubara fre- 
quency, as a function of number of iterations. QMC 
statistical error of In2T,{iTrT) at each iteration is con- 
trolled by the number of QMC sweeps while the change in 
ImY,{iT:T) between iterations is controlled by the number 
of iterations. The iteration procedure is stopped when 
the change in ImYi{i'KT) from iteration to iteration is 
less than the QMC statistical fluctuation of ImT,{iTrT), 
SQMcImT. ~ 5ImG/{ImGf. Notice that the QMC 
fluctuation of insulating solution is five to twenty times 
larger than QMC fluctuation of metallic solution as in 
Fig. |l] mostly because ImGms ^ ImGmet- This con- 
vergence criterion is more strict than the conventional 
one, which requires to stop iteration procedure when the 
norm of ||/mE(*+^^ — /ml](*)|| is smaller than a toler- 
ance. For the conventional criterion underestimates the 
contribution of the most fluctuating component to the 
convergence by taking average of self energy differences 
between iterations over whole frequency range. The con- 
vergence criteria are illustrated in Fig. |l|. 

With this convergence criterion we performed the 
DMFT-QMC calculation to find solutions of a half-filled 
single-band Hubbard model at temperatures T~l/57, 
1/64, 1/100, which turn out to be well below the crit- 
ical temperature Tc. A careful selection of initial seeding 
to feed DMFT iteration is an essential part of finding of 
the coexistence at finite temperature. At a given temper- 
ature, the best educated seeding is used to start DMFT 
iteration at initial value of interaction strength U: small 
interaction strength U for metallic solution and large U 
for insulating solution. After DMFT-QMC run we ob- 
tain convergent metallic and insulating solutions within 
a few iterations. We increase interaction strength U by 
small increment(O.Ol) for metallic solution and decrease 
U for insulating one. To reduce computational time, the 
convergent solution at the previous value of interaction 
strength U is used as initial seeding for DMFT-QMC run 
at the next value of interaction strength U. We observe 
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FIG. 1. Convergence criterion. The error bars are the 
QMC fluctuations at individual iterations, (a)lnsulating so- 
lution, (b) Metallic solution. 



convergence pattern of solutions and repeat the proce- 
dures until critical slowing down is noticeably developed. 
By plotting inverse number of required iterations for the 
convergence of insulating solution in Fig. 0(a) and of 
metallic solution in Fig. |^(b), we identify upper bound 
of UciiT) and lower bound of Uc2{T). To find lower 
bound of Uci{T) and upper bound of Uc2(T) we approach 
coexistence boundaries, UciiT) and Uc2{T), from out- 
side the coexistence region. Outside the coexistence re- 
gion only one solution exists. For interaction strength 
U{T) < UciiT) DMFT iteration which starts with in- 
sulating seeding would produce convergent metallic solu- 
tion. As U gets closer to Ud from outside the coexistence 
region, convergence rate become slower. Such critical 
slowing down is clearly captured with open symbols in 
Fig. 0(a). By approaching coexistence boundaries from 
both inside and outside the coexistence region we obtain 
upper and lower bounds of Ud and Uc2 at several tem- 
peratures: for P = 57, 2.36 < C/ci(r) < 2.4 and 2.42 < 
Uc2{T) < 2.53; for (3 = 64, 2.36 < Uci{T) < 2.41 and 
2.48 < Uc2{T) < 2.52; for /? = 100, 2.35 < Uci{T) < 2 Ah 
and 2.58 < Uc2{T) < 2.62. 

At /3 = 57 the coexistence region is narrow enough so 
as to bring the critical slowing down into play over the 
whole coexistence region. Thus at temperatures above 
/3 = 57 the coexistence region could exist as well but will 
be inapproachable numerically for it is hard to achieve 
convergent solutions. 

All known errors in DMFT-QMC should be reduced as 
much as possible in order to obtain irrefutable evidence 
of the coexistence. Those errors are the QMC statistical 
errors and the systematic errors. The latter originate 
from two sources: the finite sampling size errors in both 
direct FT and iFT and the Trotter finite size errors. 

To reduce the finite sampling size errors in direct FT 
we use cubic spline interpolation between grid points of 
Green's functions. The main idea of cubic spline interp- 
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FIG. 2. The critical slowing down at boundaries of the 
coexistence region, Uci{T) and Uc2(T). Solid(or open) sym- 
bols are obtained by approaching Uci and Uc2 from inside(or 
outside) the coexistence region along solid(or open) arrows. 



olation is to get an interpolating function that has contin- 
uous first and second derivatives at grid points. Bound- 
ary conditions are necessary to close the set of liner equa- 
tions to determine interpolation function. The usual way, 
when boundary conditions are unknown, is to put the sec- 
ond derivatives at end points to zero (natural spline) . In 
the system we treat we can obtain boundary conditions 
from Green's function. To do so, we expand Green's func- 
tion G{iuj) via k-th derivatives of G{t) at end points: 
GM=Er=o(-l)'+'(G('=nr) + G«(0+))/(*c.„)^+i. 
One can also obtain correct high frequency behavior of 
Green's functions using the moments of the spectral den- 
sity [Q. Making correspondence between these two ex- 
pansions we get the following formula for the second 
derivatives of Green's functions in the case of half filled 
Hubbard model within DMFT: G(2)(0+) -I- G^^\(3-) = 
-(C/V4 + t^) for local Green's function and G';?\o+) + 

Go (/? ) = —t^ for Weiss function. The symmetry prop- 
erty of Green's function at half filling is taken into ac- 
count: gW(o+) = i-i)''G<-''Hp-^). 

Despite of a lot of frequency points used to define 
Green's function G{iLu), the finite frequency cutoff in iFT 
results in removal of the discontinuity of G{t) at end 
points. To correct this situation wc should extract l/iu 
tail from Green's function before making iFT. To get 
correct iFT of G{iLu) one should make numerical iFT 
for G{iuj) ~ l/iuj function and analytical iFT for the 
tail 1/iuj. In the case of half-filling the correct iFT is: 
G(r) = TE„e— "-[G(zc.„) --!;]- 1/2 for r > 0. 

The corrections in direct FT and iFT are implemented 
in DMFT-QMC: correct boundary conditions in spline 
interpolation and tail correction in iFT . All artificial os- 
cillations at both linTj{iijj) and G{t) presented in DMFT- 
QMC without corrections disappears in DMFT-QMC 
with corrections. 

Having implemented the corrections in DMFT-QMC, 



we perform DMFT-QMC run to get the coexistence at 
numerically most stable points, i.e., at mid-points be- 
tween the coexistence boundaries. The statistical errors 
are reduced by increase of QMC sweeps per site upto 
3 X 10^. Several DMFT-QMC calculations are performed 
with different values oi 5t{— ^). Trotter errors are re- 
moved by extrapolation of data to the limit of 5t — > 0. In 
Fig. it is shown convergence of local Green's function 
G{t) in the limit of 5t -^ 0. Insulating solutions are all 
within the QMC fluctuation (~ 10""'). However metallic 
solutions show monotonous increases with decreasing of 
5t. In Fig. ^ it is shown convergence of double occupancy, 
{d) = {nin-^}, as a function of 6t'^. The convergence of 
metallic and insulating double occupancies in the limit 
oi St ^ confirms the coexistence within DMFT-QMC. 
In the limit oi St ^ we obtain (dmet) — 0.039 and 
(dins) ^ 0.026 for /3 = 64 and C/ = 2.44; {d„,et) ^ 0.036 
and (dins) ^ 0.023 for /? = 100 and U = 2.53. With re- 
ducing all other finite size errors except ones originating 
from Trotter approximation, the effect of Trotter errors 
on the convergence of solution is found to be small com- 
pared to errors in direct FT and iFT. 

We now use our QMC results together with previous 
numerical results obtained with other methods to draw 
the phase diagram in Fig. ||. To supplement our QMC 
data, we use the estimates of the Mott endpoint within 
QMC [p|. We also include estimates for the coexistence 
boundaries |ll| , p2[ from Exact Diagonalization(ED). At 
zero temperature, the value of Uc2 was obtained with 
great accuracy in ref M] using the projective self con- 
sistent method. Other exact methods such as the nu- 
merical renormalization group method(NRG) [H give a 
very similar value for this quantity. The region near the 
Mott endpoint cannot be studied with numerically exact 
methods for the reasons discussed in this paper. However 
we can use the Landau theory to scale the results of the 
IPT near the critical region as described in ref |^. This 
approach was used to determine the Uci (T) and Uc2 {T) 




FIG. 3. Convergence of local Green's function G{t) in the 
limit of St{— ^) — > Q. Upper(or lower) set of lines correspond 
to insulating(or metallic) solutions. 



0.05 



A 
T3 
V 

>, 

O 

c 

CO 
Q. 

O 
O 

o 

Oi 
J2 

3 
O 
Q 



0.04 



0.03 



0.02 



0.01 




OP=64 ;U=2.44 
O |3=100 ; U=2.53 



0.1 



0.2 



0.3 



5t^ 



FIG. 4. Convergence of double occupancy as a function of 
St^. Open(or solid) symbols with QMC error bars correspond 
to nietallic(or insulating) double occupancies. 



lines near the Mott endpoint as well as the location of 
the first order MIT line Uc (T) . Two crossover lines above 
the Mott endpoint |1J] were scaled in a similar manner. 
The physical meaning of the crossovers were discussed in 
ref [p5[ in connection with experimental observations in 
V2O3. Landau MIT line and crossover lines are mapped 
onto QMC phase diagram by shifting of Uc{Tc) and using 
the temperature scaling discussed in ref |Q]. 

Near zero temperature, we determine the location of 
the finite temperature MIT line, by following ref M] 
and assuming an entropy difference of ln(2) between the 
metallic and the insulatining phase. The MIT lines ob- 
tained from those two methods agree well with MIT line 
obtained by ED method fM] which we also included in 
Fig. 1^. Given the systematic finite size errors of the var- 
ious numerical approximations, it is clear that there is 
quantitative consensus ( at the level of ten to twenty per- 
cent) between numerical calculations on where the coex- 
istence region of fully frustrated Hubbard model in infi- 
nite dimensions lies. We note in passing, that the line of 
second order phase transitions reported in ref. g, is ac- 
tually crossover line at high temperature and the Uc2 (T) 
line at low temperature. 

In conclusion, we find that QMC can be used to es- 
tablish rigorously the existence of a coexistence region 
in the frustrated Hubbard model in infinite dimensions. 
Furthermore the numerical results are in full agreement 
with the predictions of the Landau theory of ref. ^ be- 
low Tc- Finally we notice that the methodology described 
in this paper should be useful when applying DMFT to 
other phase transitions in strongly correlated systems. 

We thank G. Kotliar for many useful and fruitful dis- 
cussions. This research was motivated by the stimulating 
discussions at the workshop on Theoretical Methods for 
Strongly Correlated Fermions at the CRM in Montreal. 
We have used the supercomputer at the Center for Ad- 
vanced Information Processing in Rutgers. 
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FIG. 5. Phase diagram of the paramagnetic metal to in- 
sulator transition. Coexistence boundaries (circles) with their 
uncertainties(error bars) within QMC. MIT line (solid) from 
the Landau analysis of IPX data. MIT line(dot-dashed) from 
low temperature considerations(see text). The following data 
are from other numerical works(see text for references). Sec- 
ond order critical point (square) within QMC. Uci and f/c2 
points(diamonds) and MIT line(thick long-dashed) within ED 
method. T=0 MIT point from projective method(triangle) 
and from NRG(star). Two crossovers(long-dashed line and 
shaded area) above Tc- Coexistence boundaries from QMC 
and ED methods are connected for guides to eye (dotted line). 
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